/*
    Copyright 2006-2008 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file convertstellarmodel.cc
/// Stellar model file format conversion code. Takes an output file
/// from the Leitherer et al.(1999) code and saves it as a FITS file
/// containing the stellar model parameters as well as the actual SED.
//
// $Id$

#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <cmath>
#include "counter.h"
#include "preferences.h"
#include "CCfits/CCfits"
#include "misc.h"
#include <time.h> 
#include <unistd.h> 
#include <cassert>

using namespace std;
using namespace CCfits;

// insert the cvs revision for use in the history field
const string version = "$Revision$";

void write_history (HDU & output)
{
  ostringstream history("");
  //cout << history.str();
  //history.str("");
  history << "This file was created by convertstellarmodel.cc version "
	  << version.substr(11, version.size() -13)
	  << " on ";
  time_t now = time (0) ;
  char*now_string = ctime (& now);
  now_string [24] = 0;
  history << now_string << " by user ";
  string user (getenv ("USER"));
  string host (getenv ("HOST"));
  history << user << " on host "<< host << ".";
  cout << history.str() << endl;
  output.writeHistory(history.str());
}
  
int main(int argc, char** argv)
{
  if(argc!=2) {
    cerr << "Usage:\n convertstellarmodel <parameter file>\n";
    exit(1);
  }

  Preferences _prefs;
  
  if(!_prefs.readfile(argv [1])) {
    cerr << "Error opening parameter file: " << argv [1] << endl;
    return 0;
  }
  //_prefs.print();
  cout << "Stellar model conversion program." 
       << endl << endl;

  const string sedfile = 
    word_expand (_prefs.getValue ("sedfile", string ())) [0];
  string outfile = word_expand (_prefs.getValue ("outfile", string ())) [0];
  if (!suffix_is (outfile, ".fits"))
    outfile += ".fits";

  // Conversion factor for spectral luminosity unit in model
  double fluxunitconv = _prefs.getValue ("fluxconv", double ());
  // Convert lambdaunit in model
  const double lambdaunitconv = _prefs.getValue ("lambconv", double ());
  // and time unit
  const double timeconv  = _prefs.getValue ("timeconv", double ());
  const double massconv  = _prefs.getValue ("massconv", double ());

  const string powerunit = _prefs.getValue ("powerunit", string ());
  const string timeunit = _prefs.getValue ("timeunit", string ());
  const string lengthunit = _prefs.getValue ("lengthunit", string ());
  const string massunit = _prefs.getValue ("massunit", string ());

  const bool use_lines =_prefs.getValue("use_lines", bool ());
  const bool use_yields =_prefs.getValue("use_yields", bool ());
    
  // this is the stellar mass that the SED is given for, we want it to
  // be one of what ever the mass unit is.
  double refmass = _prefs.getValue ("refmass", double ())/massconv;
  _prefs.setValue ("refmass", double (1.0),
		   "[" + massunit +"] The stellar mass for the SED");

  // this is a fudge factor that changes the mass normalization (it's
  // separate from the reference mass above only for clarity),
  // typically for IMF changes on the low-mass end. (>1 means less
  // light per mass)
  double imf_fudge =_prefs.defined("imf_fudge")?
    _prefs.getValue("imf_fudge", double ()): 1.0;
  
  // true means that the file contains log (L_lambda), as the files
  // from Starburst 99 do.
  const bool logflux = _prefs.getValue ("logflux", bool ());
  if (logflux) {
    fluxunitconv = log10 (fluxunitconv);
    refmass = log10 (refmass);
    imf_fudge = log10 (imf_fudge);
  }
    
  ifstream inf(sedfile.c_str());
  if (!inf.good ()) {
    cout << "Error: can't open input file " <<sedfile <<endl;
    exit (1);
  }
  
  // First we need to eat through the file until we get to 'TIME'
  string buffer;
  getline (inf, buffer);
  while(buffer.find ("TIME") == string::npos)  {
    getline (inf, buffer);
    if(!inf) {
      cerr << "Error: Can't find initial position in input file (TIME)"
	   << endl;
      exit (1);
    }
  }
  cout << "Reading input file " <<sedfile <<endl;
  // Next we read will be the initial time.

  // Read line file up to first time
  ifstream line_file;
  if (use_lines) {
    const string line_file_name =
      word_expand (_prefs.getValue("ewidth_file", string ())) [0];
    line_file.open (line_file_name.c_str());
    if (!line_file) {
      cerr << "Error opening hydrogen recombination line file!"  << endl;
      exit (1);
    }
    getline (line_file, buffer);
    while(buffer.find ("TIME") == string::npos) 
      getline (line_file, buffer);
    cout << "Reading hydrogen recombination lines from file "
	 << line_file_name << endl;
  }

  // Read yield file up to first time
  ifstream yield_file;
  if (use_yields) {
    const string yield_file_name =
      word_expand (_prefs.getValue("yield_file", string ())) [0];
    yield_file.open (yield_file_name.c_str());
    if (!yield_file) {
      cerr << "Error opening yield line file " << yield_file_name << endl;
      exit (1);
    }
    getline (yield_file, buffer);
    while(buffer.find ("TIME") == string::npos) 
      getline (yield_file, buffer);
    cout << "Reading yields from file "
	 << yield_file_name << endl;
  }
    

  vector<float> times, lambdas;
  vector<float> sed;
  float current_time =-1, current_lambda =-1;
  float time, lambda, flux, slask; 
  bool write_lambdas = true;
  counter c(1000);

  vector<double> line_strength (4);
  vector<float> line_lambda (5 );
  line_lambda [0] = 4861.32;// H_beta
  line_lambda [1] = 6562.80;// H_alpha
  line_lambda [2] = 12818.1;// Pa_beta
  line_lambda [3] = 21655.0;// Br_gamma
  line_lambda [4] = 1e20;
  double line_time= 0;
  int current_line = 0;
  int l = 0;

  vector<double> mass_loss;
  double yield_time=0;
  double current_mass_loss=0;

  while((inf >> time >> lambda >> flux >> slask >> slask).good()) {
    c.incr();

    // See if this time has ended
    if ( time > current_time) {
      // Put new time in times
      times.push_back(time*timeconv);

      if (use_lines) {
	// read to this time in the line file as well.  Line file
	// contains luminosities (erg/s)

	// if we need to read new line, do so until we get to current time
	while (line_time < time)
	  line_file >> line_time >> slask >> line_strength [1]
		    >>slask >>slask >> line_strength [0]
		    >>slask >>slask >> line_strength [2]
		    >>slask >>slask >> line_strength [3]
		    >>slask;
	assert(line_file.good());
	//cout << "time/h_alpha " << line_time << '\t' << line_strength [0] << endl;
	current_line = 0; 
      }

      if (use_yields) {
	// read to this time in the yield file as well. Yield file
	// contains log solar masses
	while (yield_time < time)
	  yield_file >> yield_time
		     >> slask >> slask >> slask
		     >> slask >> slask >> slask
		     >> slask >> slask >> slask
		     >> slask >> slask >> slask
		     >> current_mass_loss;
	assert(yield_file.good());
	// convert from log msun to fraction of initial mass lost
	const double ml=logflux?
	  1-pow(10.,current_mass_loss-refmass-imf_fudge):
	  1-pow(10.,current_mass_loss)/(refmass*imf_fudge);
	assert(ml>=0);
	assert(ml<=1);
	mass_loss.push_back(ml);
      }
    }

    // See if we should write a new lambda to lambdas
    if (write_lambdas && (lambda > current_lambda))
      lambdas.push_back(lambda);
    else
      write_lambdas = false;

    // perform sanity check that wavelengths are consistent
    if (lambda < current_lambda)
      l = 0;
    if (lambda != lambdas [l]) {
      cerr << "Error: Inconsistent wavelengths!"  << endl;
      cerr << "\tTime: " << time << "\tLambda: " << lambda
	   << "\tshould be: "
	   << lambdas [l] << endl;
      exit (1);
    }
    ++l;
      
    // check if we have a line
    if (use_lines && ( lambda > line_lambda [current_line])) {
      if (logflux)
	// ideally, we would like to use lambda (i+1) - lambda (i-
	// 1)/2 for the line luminosity, but we don't know the next
	// lambda...
	flux = log10 (pow (double (10), static_cast <double> (flux) ) +
		      pow (double (10), line_strength [current_line])/
		      (lambda - current_lambda));
      else
	flux += line_strength [current_line]/(lambda - current_lambda);;
      current_line ++;
    }
    
    // Write this flux point to the current sed
    sed.push_back(logflux? flux+fluxunitconv-refmass-imf_fudge:
		  flux*fluxunitconv/(refmass*imf_fudge));
 
    current_time = time;
    current_lambda = lambda;
  }

  // convert wavelength units
  transform(lambdas.begin(), lambdas.end(), lambdas.begin(),
	     bind2nd(multiplies<float> (), lambdaunitconv));
	    
  cout<< "\nSED contains " << lambdas.size() << " wavelengths and "
      << times.size() << " time steps.\n";
  cout << "covering " << times [0] << " - "<< times.back()<< " years\n"
     << "and " << lambdas [0] << " - " << lambdas.back() << " m\n";
  assert(sed.size()==lambdas.size()*times.size());

  cout << "Writing output file " << outfile << endl;

  // Now we can write the output FITS file
  if (_prefs.defined("CCfits_verbose"))
    FITS::setVerboseMode(_prefs.getValue("CCfits_verbose", bool ()));

  FITS fits (outfile, Write);
  fits.pHDU().addKey("FILETYPE", string ("STELLARMODEL"),
		     "This file contains the SED of a stellar model ");

  // Write STELLARMODEL HDU
  ExtHDU* keywordHDU = fits.addTable("STELLARMODEL", 0);
  _prefs.write_fits(*keywordHDU);
  keywordHDU->addKey ("mintime", times [0],
		      "["+timeunit+"] youngest stellar population included");
  keywordHDU->addKey ("maxtime", times.back(),
		      "["+timeunit+"] oldest stellar population included");
  write_history (*keywordHDU);

  // Write SED HDU
  // L_lambda column
  Table* sedHDU = fits.addTable("SED", 0);  
  const string llc =  (logflux?  string ("log "): string ("")) + "L_lambda";
  sedHDU->addColumn(Tfloat, llc, lambdas.size(), powerunit+"/"+lengthunit);
  Column& l_l = sedHDU->column(llc );
  l_l.write(sed, times.size(), 1);

  // time column
  sedHDU->addColumn(Tdouble, "time", 1, timeunit);
  sedHDU->column("time").write(times, 1);

  // mass loss column
  sedHDU->addColumn(Tdouble, "current_mass", 1, massunit);
  sedHDU->column("current_mass").write(mass_loss, 1);

  // Write LAMBDA HDU
  Table* lambdaHDU = fits.addTable("LAMBDA", 0);
  lambdaHDU->addColumn(Tfloat, "lambda", 1, lengthunit);
  Column& c_l = lambdaHDU->column("lambda");
  c_l.write(lambdas, 1);

  // dump
  ofstream out ("model");
  for (int i = 0; i < times.size(); ++i) {
    //cout << i << '\t' << times [i] << '\n';
  for (int j = 0; j < lambdas.size(); ++j) {
    out << lambdas [j] << '\t' << sed [i*lambdas.size()+  j] << '\n';
  }
  out << "\n\n";
  }
}

